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ABSTRACT 

We present a new code for radiation transport around Kerr black holes, including arbitrary emission and absorption 
mechanisms, as well as electron scattering and polarization. The code is particularly useful for analyzing accretion 
flows made up of optically thick disks and optically thin coronae. We give a detailed description of the methods 
employed in the code and also present results from a number of numerical tests to assess its accuracy and 
convergence. 
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1. INTRODUCTION 

Since the beginning of X-ray astronomy over 50 yr ago, there 
has been steadily growing interest in relativistic radiation trans- 
port. Because of the high energies of both photons and electrons 
relevant to these astrophysical sources, special relativistic ef- 
fects must be included in most particle interactions. Because 
the central engines of so many X-ray sources are compact ob- 
jects such as neutron stars and black holes, full general relativity 
must also be incorporated into any physically realistic code. 

Here, we present in detail for the first time the fully rel- 
ativistic Monte Carlo radiation transport code Pandurata. 3 
While this is the first formal description of the code in the 
literature, it has been under development for many years and 
has already been used in numerous publications (Schnittman & 
Krolik 2009, 2010; Noble et al. 2011; Schnittman et al. 2013). 
Pandurata shares features with numerous existing codes in the 
literature, but we believe its combination of full general relativ- 
ity, wide range of emission and absorption processes, scatter- 
ing, polarization, and optically thin/thick capabilities make it 
uniquely valuable in the rapidly evolving field of black hole 
astrophysics. Most recently, in Schnittman et al. (2013) we 
demonstrated how the code can be used to make a major step 
toward bridging the gap between global magnetohydrodynam- 
ics (MHD) simulations and real X-ray observations of accreting 
black holes. 

The literature of radiation transport in astrophysics is ex- 
tremely broad and includes scores of different techniques and 
applications. It would be a futile endeavor to attempt to give a 
comprehensive summary of the work here. Rather, we will sim- 
ply highlight a few recent contributions that are most relevant 
to the applications of interest, namely photon transport around 
Kerr black holes. By far the most common approach has been 
ray-tracing geodesic paths backward from a distant observer 
to the accretion flow, calculating a transfer function of some 


3 The name Pandurata comes from Coelogyne pandurata, a species of 
black orchid native to Borneo. Much of the core radiation transport is derived 
from the code Buttercup, developed by J.D. Schnittman for inertial fusion 
applications (Schnittman & Craxton 1996, 2000) at the University of 
Rochester’s Laboratory for Laser Energetics (LLE). In holding with LLE’s 
long tradition of naming codes after flowers, the name Pandurata was chosen 
to represent the joint heritage of black holes and laboratory radiation 
hydrodynamics. 


sort, and coupling this to some model for emission to gener- 
ate spectra and light curves. A few of the many examples of 
this observer-to-emitter approach include Rauch & Blandford 
(1994), Broderick & Blandford (2003, 2004), Schnittman & 
Bertschinger (2004), Dovciak et al. (2004), Schnittman et al. 
(2006), Noble et al. (2007), Dexter & Agol (2009), and Dexter 
et al. (2009). 

A smaller number of codes have been written with the emitter- 
to-observer approach, which may be more physically intuitive, 
but is almost always more computationally intensive, with the 
exception of Laor et al. (1990) and Laor (1991); Kojima (1991) 
who use uniform sampling of emission angles, these codes are 
generally Monte Carlo in nature, such as grmonty by Dolence 
et al. (2009) and the present work. As we will show below, 
particularly when electron scattering is included, the emitter- 
to-observer paradigm is almost essential for capturing the most 
relevant physics of the problem. 

Another feature that is relatively uncommon in these ray- 
tracing codes, but of increasing interest in the high energy 
community, is polarization. It is included in Agol & Krolik 
(2000), Dovciak et al. (2008), Huang et al. (2009), Shcherbakov 
& Huang (2011), Huang & Shcherbakov (2011), and Marin 
et al. (2012), although often only for vacuum transport and not 
including scattering. Disk polarization is treated by Laor et al. 
(1990), Matt et al. (1993), and Dovciak et al. (2008) for both 
Schwarzschild and Kerr black holes, but neglecting electron 
scattering. Dovciak et al. (2011) includes illumination from a 
source above the disk, while Dovciak et al. (2008) includes a 
cold plane-parallel atmosphere above the disk, geometrically 
thin with varying optical depth. A small number of ray-tracing 
codes also allow for non-standard black hole metrics as a way 
of testing general relativity. Krawczynski (2012) follows the 
emitter-to-observer paradigm for calculating polarized flux from 
a thermal disk and Psaltis & Johannsen (2012) describes an 
observer-to-emitter framework that can be applied to a large 
number of space-time tests such as timing, spectra, and imaging 
(Johannsen & Psaltis 2010a, 2010b, 2011, 2013). 

The body of literature including detailed scattering and 
polarization is generally restricted to flat spacetime and often 
only the most simple geometries (Connors & Stark 1977; 
Connors et al. 1980; Sunyaev & Titarchuk 1985; Haardt & 
Maraschi 1993; Haardt et al. 1994; Poutanen & Svensson 1996). 
Here, we attempt to synthesize the strengths of all these various 


1 


The Astrophysical Journal, 777:11 (17pp), 2013 November 1 


ScHNITTMAN & KROLIK 


codes into a single flexible radiation transport tool for analyzing 
both global MHD simulations and also simpler toy accretion 
models. The ultimate goal is to produce concrete predictions 
that can be compared directly with the large and continually 
growing body of X-ray observations of accreting black holes. 

Precisely because of the large interest in this topic, we give 
here a comprehensive description of the technical components 
of our radiation transport code Pandurata. We hope that the 
techniques outlined below will be valuable to others who are 
interested in developing similar (or even better, more powerful) 
tools. We also present the results from a suite of simple 
numerical tests to verify the code, thus lending support and 
increasing our confidence in earlier work based on Pandurata. 

2. LOCAL ORTHONORMAL FRAMES 

The most general input for Pandurata is a body of tabu- 
lated data including the extrinsic fluid variables density, tem- 
perature, magnetic field, and the four-velocity at each point in 
a three-dimensional volume. Multiple data slices in the time 
coordinate allow for studies of variability. The coordinates are 
Boyer-Lindquist for a Kerr black hole with mass M and dimen- 
sionless spin parameter a/M. The fluid variables are given in 
physical cgs units for a specific black hole mass and accretion 
rate. The source of the data is quite general and Pandurata 
has been used successfully to analyze simulation data from the 
relativistic MHD codes Harm3D (Noble et al. 2011; Schnittman 
et al. 2013) and GRMHD (Schnittman et al. 2006), as well as 
two-dimensional hydro simulations (Schnittman & Rezzolla 
2006) and analytic models for the accretion disk and corona 
(Schnittman & Krolik 2009, 2010). 

We adopt a (— + ++) metric signature and a convention where 
Greek indices run from 0 to 3 and Latin indices are restricted to 
spatial coordinates from 1 to 3. The coordinate metric is given 


by (Boyer & Lindquist 1967) 





(~ 

1 11 

-or + 

0 

0 

— com 2 \ 

8/j.v — 

0 

P 2 / A 

0 

0 


0 

0 

P 2 

0 


\ 

—com 2 

0 

0 

m 2 ) 


This allows for a relatively simple form for the inverse metric: 
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In geometrized units with G = c = 1, we have 
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Figure 1. Fluid density (top) and temperature (bottom) profiles for a slice of 
Harm3D data in the (r, z) plane, taken from the ThinHR simulation (Noble 
et al. 2010, 2011). Contours show surfaces of constant optical depth with 
r = 0.01, 0.1, and 1.0. Fiducial values for the black hole mass M = 10 Mq and 
luminosity L = 0.1 Z,Edd were used, as described in Schnittman et al. (2013). 
(A color version of this figure is available in the online journal.) 

2.1. Simulation Data 

We briefly describe here the format of data from the Harm3D 
MHD simulations. Similar data can be generated from GRMHD 
and any analytic model can be understood as a subset of the 
full tabulated simulation data. As described in greater detail in 
Noble et al. (2011) and Schnittman et al. (2013), the first step 
in post-processing the simulation data is to convert from code 
units of density and local dissipation to cgs units of density 
and temperature. Given the density everywhere, we integrate 
the optical depth along paths of constant ( r , </>) coordinates 
starting from both 0 = 0 and 9 = jt to get r top (r, 9, cp) and 
TbotO", $, </>)• The disk midplane can then be defined as the 
surface 0 n iid(r <P) where r top (r, 0 m id, <P) = r bo t (r, 0 m id, </>)• When 
r(r, 0 m id, 4>) > 1, the disk is optically thick and we define a 
top and bottom photosphere 0(r, <p) such that r top (r, ©top, 4>) = 
r b „t(r, ©hot, <P) = 1. In Figure 1, we show a slice in the 
(r, z ) plane of simulation data from the Harm3D “ThinHR” 
run (Noble et al. 2010). The local temperature is represented 
by the logarithmic color scale and the contours show surfaces 
of constant r. In Figure 2, we show a three-dimensional 
representation of the photosphere surface 0 top (r, <p) for the same 
simulation data. 

From the photosphere surfaces, thermal photons are launched 
into the optically thin corona above and below the disk. Because 
the opacity within the disk is usually dominated by electron 
scattering, the seed photons are emitted with the limb-darkening 
and polarization dependence on angle given by Chandrasekhar 
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Figure 2. Three-dimensional representation of the disk photosphere surface 0 top (r, 0), along with the local spatial tetrad definitions of Equation (18a). The simulation 
data are the same as in Figure 1. The color scale is a linear representation of the disk’s local thermal temperature. The labels (r; , 0/) correspond to coordinates of the 
computational grid boundaries. 

(A color version of this figure is available in the online journal.) 


(I960). The spectrum is that of a diluted blackbody with 
temperature T e g and hardening factor/: 


I v 


1 

T a 


B v (f T eB ), 


(4) 


where B v is the blackbody function. We take / = 1.8 (Shimura 
& Takahara 1995) and the local effective temperature is 
given by 


T e s(r, (p) = 


Hr, <j>) 


1/4 


(5) 


where 2 T(r, </>) is the total integrated flux from the optically 
thick part of the disk (the factor of 2 comes from the fact that the 
flux is emitted equally from the top and bottom photospheres). 
As shown in Figure 1 , the gas has a constant temperature inside 
the disk for a given (r, dp), due to the high level of thermalization 
caused by the large optical depth. 

Synchrotron and bremsstrahlung seed photons can also be 
generated in the coronal regions, in which case we use an 
unpolarized, isotropic distribution function for the emission 
angles, as measured in the local fluid frame and by the angle- 
averaged formulae of Mehadevan at al. (1996). Because the 
current version of Pandurata does not include angle-dependent 
synchrotron emission, all polarization comes from electron 
scattering. Future work will include polarized synchrotron seeds 
as well. Due to the high temperatures and low densities of the 
coronal regions, the net power in the coronal seed photons is 
typically much lower than that of inverse Compton scattering 
from the thermal seeds coming from the disk (Schnittman et al. 
2013). 


2.2. Photosphere Tetrads 


e” ) = S^, where S is the usual Kronecker delta. Note that 
the coordinate basis vectors are not normalized and not even 
orthogonal in the Kerr metric: 

e O) ' e (v) = gafie ( ^e (v ) = g ^ v . (6) 

Einstein’s Equivalence Principle, one of the bedrocks of 
general relativity, states that an orthonormal basis (a “tetrad”) 
can be defined at any point in space. In fact, an arbitrary number 
of tetrads can be defined at any point and are all related by 
Lorentz boosts and/or rotations. One particularly useful tetrad 
in the Kerr metric is that of the zero angular momentum observer 
(ZAMO; Bardeen et al. 1972). We denote the ZAMO frame with 
jl labels. It can be constructed from the coordinate basis by 
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Any vector can be represented by its components in different 
bases: 

u = e (lLt) « M = e (ji )/ (8) 

and the components are related by a linear transformation Et : 
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We begin with a short discussion of notation. As stressed 
in Misner et al. (1973), vectors are invariant geometric objects 
independent of coordinate system and we represent them with 
bold font u, while the components in a specific basis are 
represented with italics: u 11 . We adopt a naming convention 
such that the components of a vector in the coordinate basis 
are represented by /i and in the local fluid frame by (i. The 
basis vectors themselves are labeled with (//). For example, the 
coordinate basis is spanned by the vectors with components 
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At each point on the photosphere, we define a tetrad in the 
comoving fluid frame (designated with sub/superscripts jl ) such 
that the time coordinate is in the direction of the fluid four- 
velocity: 


et 

(?) 


( 11 ) 


photosphere two vectors tangent to the disk surface through the 
following process. We start with the coordinate-based vectors 


dr* = 


90 

0, A r, — A r, 0 
dr 


(14a) 


In our notation, this equation says that the four-vector tangent 
to the world line of an observer moving with the fluid can 
be expressed in the Boyer-Lindquist coordinate basis with 
components /i or in the local frame with components / 1 with 
= [1,0, 0, 0]. The spatial basis vectors in the fluid frame 
are constructed via a method similar to Beckwith et al. (2008), 
including a slight modification to ensure the right-handedness of 
the basis such that eg) is in the —0 direction. For completeness, 
we reproduce those definitions here: 


and 


dtp" 


90 

0, 0, — A0, A</> 

o(p 


(14b) 


where A r and A cp are the differential sizes of the fluid cell in 
question. 4 Next, we subtract off the components parallel to eg) : 


dr' = dr - (dr e (f) )e (?) (15a) 
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and 

d0' = d0 - (d</> • e (? ,)e (?) . (15b) 

When dr' and A(j> are projected onto the fluid frame, they 
will have only spatial components and will be tangent to the 
photosphere. In this basis, we can easily construct the normal 
vector by taking the three-vector cross product: 

dz k = e\,dr' 1 d/> . (16) 

l J 

This procedure has the added advantage of giving the proper 
area of the photosphere patch subtended by the vectors dr' 
and A<j>' by d A = |dz|. This formula for d A will be helpful 
for determining the amplitude of emitted flux from each patch 
of the disk, since the emission function is typically defined in 
the local fluid frame. Because dr' and A<j>' are not generally 
orthogonal, we also define the dx and dy tangent vectors by 
dx = dr', dy' — 0, and 


dy k = e\,dz'dx^ . 
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To complete the tetrad, we simply need to normalize the 
differential basis vectors. Returning to the coordinate basis, we 
have: 


(13a) 

M _ H 

e ff) ~ e (i) ~ 

(18a) 

(13b) 

eg) = dx ll /(g a pdx a dx l} ) 1/2 

(18b) 


4 ) = /(g<*pdy a dyP) 1/2 

(18c) 

(13c) 

eg) = ±dz^ l(g a pdz a dz p ) 112 . 

(18d) 


From this tetrad basis, any other tetrad in the fluid frame can 
be constructed from a simple rotation of the spatial basis vectors. 
We take as our preferred basis (now labeled with e(g)) one in 
which eg) is normal to the photosphere surface. Whether we are 
using simulation data or an analytic model for the disk surface, 
the photosphere is described by a two-dimensional surface 
on the top and bottom of the disk: 0 tO p(r, </>) and 0h„ t (r, <p). 
From these functions, we can construct at each point in the 


The ± in the definitions for e ( vl and eg) are chosen for the top 
(+) and bottom (— ) photosphere surfaces so that the spatial basis 
vectors are oriented in a right-hand fashion and to ensure that eg) 
points away from the disk body. In Figure 2, we show how these 
tetrad basis vectors are oriented on the photosphere surface. 


4 For example, the ThinHR simulation uses Ar/r = 0.004 and A (j> = tr/128. 
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2.3. Coronal Tetrads 


3. RAY-TRACING 


In addition to launching photons front the disk surface, 
we often want the option of including seeds from within the 
corona, due to thermal bremsstrahlung, cyclo-synchrotron, or 
other radiation processes. Analogous with the comoving surface 
element defined above for disk emission, for coronal emission 
we need to define a volume element and associated tetrad at each 
point in the simulation space. Like the tetrads defined above, 
the time axis is defined along the local fluid four- velocity u ,L . 
However, unlike the surface tetrads, the volume tetrads have 
no preferred orientation,' 1 so we can simply use the spatial 
coordinate vectors projected onto the space orthogonal to u !> : 


dr^ = [0, A r, 0, 0] , 

(19a) 

dd" = [0, 0, A6>, 0] , 

(19b) 

dtp" = [0, 0, 0, A 0] , 

(19c) 

dr' = dr — (dr ■ e ( /j) , 

(20a) 

d 6' = d6 - «W ■ e (?) ) , 

(20b) 

d0' = d0 - (d0 • e (f) ) . 

(20c) 


The proper volume element subtended by these vectors is given 
by the three-vector triple product in the local fluid frame. While 
there is no real preferred orientation for the spatial axes, we still 
need to go through the process of defining some orthonormal 
basis to project the above vectors and thereby calculate vector 
products. In practice, we set eg) along the dr' direction: 

dx* = dr' 11 , (21a) 

then set the y-axis roughly along the </; coordinate direction 


3.1. Geodesics 

The ray-tracing portion of Pandurata integrates the geodesic 
trajectories of photons in the Kerr metric. From the tetrad frames 
defined in the previous section, the initial direction of a photon 
is selected from an isotropic distribution in the emitting fluid 
frame (limited to a hemisphere in the case of an optically thick 
photosphere surface). 

The geodesic integrator is the same as that described in 
Schnittman & Bertschinger (2004), based on a Hamiltonian 
approach. Because the Kerr metric is stationary, the momentum 
conjugate to the time coordinate t is conserved and corresponds 
to the (negative) specific energy of a particle ( m 2 = 1) or 
photon ( m 2 = 0). We can replace the affine parameter with 
the coordinate time and write the Hamiltonian as 


H(x l , pi) = -p 0 


will] equations of motion 


8? Pi , 

g ,J PiPj+m 5 2 

1 

<N 

O 

g 00 
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dx‘ 3 H\ 

dt dpi 


(24a) 


dp; _ dHi 
dt dx l 


(24b) 


In Boyer-Lindquist coordinates, the Hamiltonian can be 
written thus: 


H(r, 0, <p, p r , p g , Pcp) 
= cop(/) + a 


A 2 , 1 2 , 1 2 , 2 

-p Pr + -pPe+ —p<t, + m 
p- P A UJ A v 


1/2 


(25) 


dy ic = e\,dx ri d0 ] (21b) 

and the z-axis normal to both: 

dz^ = z\jdx 1 dy’ . (21c) 

As above for the photosphere tetrads, the final step is to 
normalize all the basis vectors: 


= dx^ /(g a pdx a dx^) 1/2 (22b) 

e$) = dy 1 * / (gapdy 01 dyP) 1/2 (22c) 

eg) = dzy(gapdz a dz p ) l/2 . (22d) 

Unlike the photosphere case, since there is no “top” or “bottom” 
in the corona, we need not be concerned about the orientation 
of the eg) vector and simply require a right-handed (x,y,z) 
convention. 

5 For some specialized emission models, such as optically thin synchrotron, it 

may be convenient to choose a special orientation, e.g., with the eg) basis 

rotated to lie along the local magnetic field vector. 


using the same notation defined above in Equations (3a)-(3e). 
Because the metric and thus the Hamiltonian, is axisymmet- 
ric, p,i, is also an integral of the motion. We are thus left 
with five coupled first-order ordinary differential equations for 
(r. 6, (j> . p r . p e ). The third integral of motion, Carter’s constant 
(Carter 1968) is 

Q = Pe+ cos 2 0 [ a 2 (m 2 - pi) + pj/ sin 2 0] (26) 

and is used as an independent check of the accuracy of the 
numerical integration. 

For the numerical integration of geodesics, we use a fifth- 
order Cash-Karp algorithm with an adaptive step size (Press 
et al. 1992). In Figure 3, we show the accuracy of the integrator 
by plotting the average deviation in the Carter constant for a 
selection of photons around a black hole with a/M = 0.99, as 
a function of step segments. We typically set the tolerance at 
10~ 8 per step, which we find allows sufficient sampling of the 
fluid variables near the black hole. Because of the frequent table 
look-ups required when using simulation data, there is little to 
be gained by using more advanced integration techniques such 
as Bulirsch-Stoer or the semi-analytic approaches of Rauch & 
Blandford (1994) or Dexter & Agol (2009) that calculate the 
geodesic endpoint in a single integral evaluation and are more 
appropriate for vacuum transport. 
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N steps 

Figure 3. Convergence of the geodesic integrator, as determined by the accuracy 
of the conserved quantity Q. As expected, we find fifth-order convergence for 
the Cash-Karp adaptive time step integrator. 

3.2. Polarization 


3.3. Photon Packets 

Because the geodesic photon trajectories are independent of 
photon energy, we can significantly improve the efficiency of the 
Monte Carlo calculation by tracking large numbers of photons 
simultaneously covering a range of energies. We call these 
computational entities “photon packets,” which are analogous 
to the “superphotons” of Dolence et al. (2009), except for the 
fact that theirs are monoenergetic and ours are broadband. We 
also assign a single polarization angle and degree to the entire 
photon packet. This is an approximation that works well for 
vacuum transport and coherent scattering, but will break down 
when including scattering at high energies hv > m e c 2 as the 
electron cross section becomes more energy-dependent. 

Each photon packet is weighted by a number of geometric 
emission factors. For example, a photon packet emitted from a 
small patch of optically thick, scattering-dominated accretion 
disk would have a spectrum of 

F 7 = 4 B v {f T e ff) ± f limb (& cm ) cos 9 em dA dQ. , (30) 

/ u 


Pandurata is also capable of polarized transport along 
geodesics. The polarization vector is a space-like vector normal 
to the photon direction. For a photon with wavevector k. the 
polarization vector f is constrained by f • f = 1 and f • k = 0 
(Connors et al. 1980). The vector f is parallel transported 
along the geodesic path: V^f = 0, but instead of explicitly 
solving the parallel transport equation, we can take advantage 
of the complex-valued Walker-Penrose constant /c wp (Walker & 
Penrose 1970; Connors & Stark 1977). 

After solving for the wavevector k 11 along the geodesic path, 
/c wp is given at any point by 

Kwp = [(k‘f r - k'-f ) + a sin 2 9{k r f+ - k* f) 

- i[(r 2 +a 2 )(k<’f e - f*k e ) - a(k'f - k e f)} sin0] 
x (r — ia cos 9) . (27) 

Combined with the normalization factors f • f = 1 and f • k = 0, 
we have four linear equations for the four components of f . 
Because k is a null vector, we can always redefine f by a multiple 
of k: f' = f + /.k and thus write the polarization vector as 

/** = [0, cos \j/e\ + sin i/re^] (28) 

for some space-like basis vectors ei and e2 normal to k. 

The degree of polarization 5 ^ 1 is invariant along the ray 
path. When interacting with a distant detector or scattering off 
an electron in the fluid frame, it is convenient to employ the 
classical Stokes parameters /, Q, and U. In the (ei, 62) basis, we 
can write 


X = Q/ 1 = 8 cos 2i[r , 

(29a) 

Y = U / / = 8 sin2i jr . 

(29b) 


One of the main advantages of this approach is that the Stokes 
parameters for each photon can simply be added at the detector, 
quite useful in a Monte Carlo calculation. Furthermore, I(v), 
Q(v), and U(y) can all be written in units of spectral density, 
which is the standard observable for many real detectors. 

For photons emitted at an angle 9 em to the normal of a 
scattering-dominated surface, we use the results of Chan- 
drasekhar (1960) to obtain the initial polarization amplitude 
(ranging from <5(0 em = 0°) = 0 up to 8(9 em = 90°) = 0.12) and 
direction (parallel to the disk surface). 


where F v is a function that has units of spectral luminos- 
ity [erg s -1 Hz' 1 ]. Here, / is the same hardening function 
introduced above in Equation (4), cos 9 em is a geometric fac- 
tor for emission from an optically thick surface, is a limb- 
darkening function given by Chandrasekhar (1960), d A is the 
proper area of the emission region (see Equation (16) above), 
and di 2 = 2jr//V ph is the proper solid angle of a hemisphere 
sampled evenly by A p h photon packets. Lastly, 1 /m' = dr /dt is 
a relativistic correction factor to convert from time in the emis- 
sion frame to that in the coordinate or distant observer frame. 

In order to account for the spectral redshift, we store both F 
and v at a set of discrete points. When transforming from the 
emitter to observer frames, F is invariant (units of s _1 and Hz -1 
cancel), 6 while v transforms as follows. If the photon packet is 
emitted in a frame with fluid four- velocity u lx (e m) and photon 
four-momentum l M (em) and observed in a frame with m 11 fobs) 
and C/obs), then we can write the redshifted frequencies as 


v(obs) = v(em) 


M l ’(obs)A:„(obs) 

M^(em)k lU (em) 


(31) 


Whenever the photon packet scatters off the disk or an electron in 
the corona, the frequencies v* are updated and the old “observed” 
frame becomes the new “emitted” frame. When the photon 
packet reaches an observer at infinity, iT/obs) = [1,0, 0,0] 
and the well-known redshift relation is reproduced. 

For this distant observer, the angle of polarization 1 [r is 
measured by projecting f onto the (ei , e2) = (e^, —eg) basis. For 
an observer oriented with the black hole spin axis projected in the 
vertical direction, ifr = 0 corresponds to horizontal polarization 
(Schnittman & Krolik 2009, 2010). Given \l/, S, and F v , the 
spectral luminosity form of the Stokes parameters are simply 


Q v = F v 8 cos 2\[r , (32a) 


U v = F v 8 sin2i/r , 


(32b) 


6 For a discrete function / ;. the number of photons emitted per 
coordinate-frame second between v,- and v, + dvj is F, dvi/{hvi), where h is 
Planck’s constant and v,- are measured in the local emission frame. Because v,- 
and dvi transform the same under Lorentz transformations, F, is invariant. 
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where Q and U are related to the Stokes parameters Q and U 
by a factor of ( F/I ). After summing over a large number of 
photon packets, we then invert Equation (32) and return to the 
<5(v), i/r(t>) representation. 

3.4. Emission and Absorption 

Along each geodesic path, we can also include local emis- 
sion and absorption processes such as bremsstrahlung or syn- 
chrotron. This is the predominant method for generating light 
curves and spectra in codes that shoot photons backward from 
a distant observer (Broderick & Blandford 2004; Schnittman & 
Bertschinger 2004; Schnittman et al. 2006; Noble et al. 2007, 
2009; Dexter & Agol 2009). In the fluid frame, the radiation 
transport equation is given by 

dl v 

— = jv - ot v I v , (33) 

ds 

where ds is the differential path length and j v , and a v are the 
spectral intensity, emissivity, and absorption coefficient of the 
local fluid, respectively. The absorption coefficient is related to 
the opacity k v through the density p: a v = pK v . Defining the 
optical depth r v through 

d r v = a v ds, (34) 


4.1. Coronal Scattering 

The first step in the scattering process is to determine whether 
a scattering event takes place at all. To do this, we transform 
into a local inertial “lab” frame, generally taken to be the ZAMO 
frame discussed above in Section 2.2. In this frame, the photon 
moves a distance of dl 2 = ryjdx 1 dx J in a single geodesic 
integration step dt. Then, the total optical depth to scattering 
along the path is 

dt — dl K es Pl a b — dl /C es Pfluid — ? (37) 

Hab 

where the last equality comes from the invariance of va v 
(Rybicki & Lightman 2004), with the absorption coefficient 
a v = K es p for electron scattering opacity. Given dr (typically 
much less than unity), the probability of scattering is P = 

1 - e~ dz . 

When a photon does scatter off a free electron, we carry 
out the scattering calculation in the electron’s rest frame. This 
requires two coordinate transformations: from the coordinate 
basis (denoted with /i super/subscripts) to a fluid-frame tetrad 
(/I), and then a Lorentz boost from the fluid frame to the 
electron’s rest frame ifi'). The transformation from coordinate 
basis to corona fluid frame is the same as given above in 
Section 2.3. In the fluid frame, the electron velocity ft — v/c is 
taken from an isotropic Maxwell-Juttner distribution 


the transfer equation can be written as 
dl v 

-j^ = S v - I v , (35) 

Cl Ty 

where the source function is defined as S v = j v /a v . 

Both I v and .S',, have the same properties under Lorentz 
transformations, namely I v /v 3 and S v /v 3 are both invariant. 
Other invariant terms are the optical depth r„, va v . and j v /v 2 
(Rybicki & Lightman 2004). Thus, if we can solve the non- 
relativistic radiative transfer equation (Equation (33)) in the 
local fluid frame, then in any other inertial frame (e.g., the 
ZAMO tetrad), the special relativistic version can be written 



Here, the fluid frame (where j v and a v are defined) is the primed 
frame and the “lab frame” is unprimed. 

The above analysis, while quite useful for special relativistic 
flows in the locally flat ZAMO basis, ignores all general 
relativistic effects of curved spacetime around the black hole. 
To include these effects, we need only to shift the frequencies v, 
from one geodesic step to the next, due solely to the gravitational 
redshift, and we can treat each computational step as a new 
observer relative to the previous step, as in Equation (31). 

4. SCATTERING 

We allow for two types of scattering in Pandurata: Compton 
scattering off free electrons in the corona and scattering off 
an optically thick disk (which in turn is characterized by 
repeated scatterings in the relatively cool atmosphere). Because 
electron scattering conserves photon number, our photon packet 
approach is ideal for modeling these processes. 


f(v) = 


y 2 P 

9t K2(1/9t) 


ex P (~Y/0t) , 


(38) 


where y = 1 /sj 1 — ft 2 , 6j = kT /m e c 2 , and Ko is the modified 
Bessel function. We refer the reader to Appendix B for a 
description of our algorithm for generating a Monte Carlo 
sample of velocities that satisfy Equation (38). 

Following Misner et al. (1973), we construct a generic 
Lorentz boost in the direction of the electron four-velocity 
uk- — [y, yfinl]: 


u* L = [y,yPn 1 ] (|n| = 1), 

A',- = Y, 

A A = A J f = —fiyn J , 

K r = aK = (y - 1 )n j n k + #. (39) 

k J 


The photon momentum in the electron frame is thus given by 
pC = A *pk. 

Without loss of generality, we can carry out one more 
transformation and define the initial photon direction to lie along 
the z-axis in the electron frame. The x—y plane is decomposed 
into e j and S 2 , where the initial polarization is aligned with e\. 
The scattered radiation k f makes an angle 0 with e \ and 0 with 
k,. as shown in Figure 4. For unpolarized incident light, we can 
define e\ to lie in the plane of k, and k/, with 0 + 6 = 90°. 

For photons polarized along e \ , the angle-dependent cross 
section a{6) is given by the dipole scattering formula (Rybicki 
& Lightman 2004): 


da 

dQ. 


= Tq sin 2 0 = ^(sin 2 0 + cos 2 6 cos 2 0) , 

pol 


(40) 
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Figure 4. Schematic of the scattering geometry in the electron frame. The 
incoming radiation is polarized along the direction. The scattered radiation 
k f makes an angle 0 with e\ and 9 with k, . When projected onto the e\ — 
plane, k / makes an angle 0 with ei . 



Figure 5. Definitions of polarization axes in pre- and post-scattering coordinates. 
k| , ky , g|| , and ej, are all in the same plane, while ej_ and e' ± are normal to that 
plane. 


where </> is the standard azimuthal angle measured with respect 
to e i . Here, the classical electron radius r () is given by 

e 2 

r 0 = = 2.82 x 10 “ 13 cm. (41) 

m e c 2 

For photons scattering in the k, -£2 plane, the cross section is 
constant: da { 0 = 7t/2)/dQ. — rf } . For unpolarized light, we 
define e | as lying in the scattering plane, so the angle with 
respect to ei is n/2. Because unpolarized light is an equal 
combination of ei - and S 2 -polarized photons, we can reproduce 
the familiar cross section for unpolarized scattering: 



unpol 


l di 2 


1 

2 

^ r o(l + cos 2 9) . 


pol 


/ da {n/2) 

V dU. 


pol 


(42) 


For an arbitrary polarization degree <5, the cross section can 
be written as the sum of unpolarized light with weight (1 — 5) 
and purely polarized light with weight 5 : 


— = -rn( 1 — 5)(1 + cos 2 9) + /-q 5(1 — sin 2 9 cos 2 (p) 
d Q 2 

= —rg [(1 + cos 2 9) — 8 sin 2 9 cos 2 <p] . (43) 


Given the angle-dependent cross section, we can either pick 
the scattering angles {9, (p) directly from a distribution function 
derived from Equation (43) or, alternatively, we can pick the 
angles from a uniform distribution and give the scattered flux 
a weight based on the cross section. We compare these two 
methods in Appendix C. 

Once the new photon direction is determined, we need to 
calculate the angle and degree of the post- scattered polarization. 
Here, we follow the Rayleigh matrix method described in 
Chandrasekhar (1960). We define yet another coordinate system 
with £3 parallel to k ( , ey in the scattering plane defined by 
k, and ky and e j_ normal to that plane. Likewise, we define 


a post-scatter frame with e ' 3 parallel to k/, e' ± — e±, and 
g| in the scattering plane, but normal to k/ (see Figure 5). 
In this frame, the initial polarization vector can be written 
f, = cos i//£|| + sin \l/e L and the final polarization is fy = 
cos i/r'ey + sin i p'e' ± . 

The standard Stokes parameters are given by the intensity I, 
Q = SI cos 2i jr, U — 81 sin 2 \p, and V — 0 (electron scattering 
never leads to circularly polarized light). We further define 

/,! = l(/+®= 1(1 -5)7 + 5/ cos 2 ^ (44a) 

!± = \{I ~ Q)= ~ S)I + 81 sin 2 f (44b) 

l=[I h I±,U,V] (44c) 

and the Rayleigh scattering phase matrix 


/cos 2 9 

0 

0 

0 A 

0 

1 

0 

0 

0 

0 

cos 9 

0 

V o 

0 

0 

cos 9/ 


Then, the scattered Stokes parameters are given simply by 
I' = RI, /' = 7|j + l' ± , and Q' = /jj — l\ . Note that the cross 
section (Equation (43)) can be reproduced by writing 

I' = cos 2 <9/|| + I± = -(1 — 5)7(1 + cos 2 9) 

+ 57(cos 2 9 cos 2 xp + sin 2 ip) , (46) 


giving 

7' 1 

— = -(1 — 5)(1 + cos 2 9) + 5(1 - sin 2 9 cos 2 ip) , (47) 

now with ip taking the place of (p from Equation (43). 
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Lastly, ff is constructed by 

_ s/Q ' 2 + U ' 2 


(48a) 


ijr 1 


1 

2 


tan ~\U\ Q ') 


(48b) 


if = cosi/r'e|| +sinx[r , e l ± . (48c) 

At this point, the polarization vector and photon four-momentum 
are transformed back into corona fluid frame, then to the 
coordinate frame, and then the geodesic propagation continues 
as before, until the photon packet scatters again, is absorbed by 
the black hole, or reaches a distant observer. 

During this scattering process, the photon packet’s array 
of frequencies had to be adjusted three times: once when 
transforming from the fluid frame to the electron rest frame, 
once when losing energy to the electron recoil, and once 
when transforming back to the fluid frame. The first and last 
transformations are simple Lorentz boosts and the frequency 
scales like the photon energy: v'/v = p' ///, with p' = A'^p 1 *. 
For the scattering losses, we need to scale the frequency bins 
such that the number of photons in each bin is conserved, while 
losing energy according to the Compton recoil formula: 


E f = 


Ej_ 

1 + A(1 - 

m e c z v 


cos 9) 


Thus, the frequency scales like 

v' 
v 

and the size of each bin scales like 

dv' 
dv 


hv 

1 + -(1 — COS0) 

m e c 2 


1 + 


hv 


m e c ^ 


r ( 1 - COS0) 


(49) 


(50) 


(51) 


The number of photons in each bin dN v = F v dv/(hv ) is 
conserved in the scattering event, so we find that the effect 
of Compton recoil on the spectral luminosity is 


F' 

V 

F v 



hv 

-(1 cos 6 s ) . 

m e c 


(52) 


At very high energies hv m e c 2 , this leads to a “pile up” 
of photons and large peaks in the photon packet spectrum. 
In reality, this effect would be mitigated by incorporating 
Klein-Nishina cross sections, which decrease with energy, yet 
are incompatible with our photon packet approach that treats all 
photons as identical regardless of frequency. 7 In practice, we are 
generally interested in problems where the characteristic photon 
energies are significantly below m e c 2 , so the photon pile up is 
rarely an issue. 

While some energy is lost to Compton recoil in the electron 
frame, the more typical effect is inverse Compton scattering, 
where energy is transferred from the electrons to the photons. 
For electrons with Lorentz factors y in the fluid frame and 

7 Relativistic corrections to the cross section are on the order of 10% at 
hv = 30 keV and 50% at 350 keV. Our spectra based on Thomson cross 
sections should be at least this accurate in the corresponding energy range. 


low-energy photons with hv/m e c 2 <<C y 2 — 1, the ratio of 
energies of the photons before scattering, in the rest frame of 
the electron, and after scattering is roughly 1 : y : y 2 (Rybicki 
& Lightman 2004). For coronal electrons with temperature 
~140 keV, low-energy seeds will, on average, double their 
energy after every scattering event, making inverse Compton 
scattering a very efficient radiative process. 

4.2. Disk Scattering 

At each step along the geodesic trajectory, we determine 
whether or not the photon packet has crossed the photosphere 
surfaces 0 lIip (r, </>) or 0 ho t(r, <•/;). If it has crossed this boundary, 
we follow a procedure similar to that described above for coronal 
scattering. First, we use the conserved quantities K wp , f ■ f = 1 
and f • k = 0 to solve for the polarization vector f in the 
coordinate frame. Then, we transform f and k into the local 
fluid frame of the photosphere tetrad e^), with e,j) normal to the 
disk surface, as in Equation (18a). 

In this frame, the scattering off the disk surface is calculated 
using the analytic expressions for reflection off a diffuse, 
semi-infinite plane, derived by Chandrasekhar and given in 
Table XXV of Chandrasekhar (1960). As in Equation (44) 
above, we can write the incoming photon beam as a vector 
of Stokes parameters for the flux I = [/y, l±, U ] ( V =0 for 
linearly-polarized light, the only relevant case for our scattering- 
dominated systems). Then, the outgoing intensity is given by 

//"\ 

QS(/x, cp; /x 0 , (po) I Cl j , (53) 


I'O t,tp) 



whore ( /i (h (p 0 ) are the incident angles in the fluid frame and /x o = 
\ki\, (/x, (p) are the outgoing angles and Q and S are the transfer 
matrices defined in Section 70.3 of Chandrasekhar (1960). 
Unlike the coronal scattering case, where we use the differential 
cross section (Equation (43)) to determine the post-scatter 
angles, for diffuse reflection off the disk, we simply choose 
a random angle (/x,<p) from a uniform distribution and then 
weight the outgoing intensity by /'// from Equation (53). Thus, 
any individual reflection does not conserve photon number, 
but the angle-averaged process does. From I', we are able to 
reproduce 8', if/' , and thus f' and k' as above, which are then 
transformed back into the coordinate frame and continue their 
geodesic propagation through the corona. 

This method for diffuse reflection can be checked against 
coronal scattering experiments where we scatter incoming pho- 
tons off a semi-infinite plane of free electrons. We find excellent 
agreement between the analytic and numerical approaches, as 
shown below in Section 5. 

As with the coronal scattering, high-energy photons can lose 
energy to Compton recoil off the electrons in the cool disk, 
leading to the reflection hump seen in many active galactic 
nucleus (AGN) observations (we are able to reproduce this fea- 
ture above ~30 keV and compare with similar results calculated 
in George & Fabian 1991, who calculated the reflection spec- 
trum off a cold disk when irradiated with an external power-law 
flux). While this process is technically angle-dependent, as a 
simplification, we average over all incoming and outgoing an- 
gles, as well as the number of individual scatterings typically 
responsible for diffuse reflection ( /V sca t ~ 3 in the Thomson 
regime) and use the recoil formula 


v' 

V 



hv 
m e c 2 


-l 


(54) 
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Figure 6. Comparison of the observer-to-emitter and emitter-to-observer ray- 
tracing paradigms for a relativistically broadened emission line. The disk 
extends from if out = 15 M all the way to the horizon and the emissivity profile 
scales as r -3 . 

(A color version of this figure is available in the online journal.! 


This energy lost by the photons can then be reprocessed by 
the disk and emitted at thermal energies. In future work, 
when including more accurate Klein-Nishina cross sections 
that decrease with energy, we expect that high-energy photons 
incident on the disk will penetrate deeper into the disk, then 
lose energy to Compton recoil before having to scatter back out 
of the disk atmosphere. Thus, we expect /V scat will increase and 
thus the energy losses would be even greater. 

5. NUMERICAL TESTS 

Here, we present a number of test problems to verify 
Pandurata’s accuracy and reliability. We begin with vacuum 
transport of geodesics from the disk to a distant observer. To test 
the tetrad construction methods outlined in Section 2.2, we cal- 
culate the relativistic broadening of iron lines from a thin disk 
around a Kerr black hole, comparing the emitter-to-observer 
and observer-to-emitter paradigms. The observer-to-emitter ap- 
proach is well-known in the literature (Rauch & Blandford 1 994; 
Broderick & Blandford 2003, 2004). It is also relatively straight- 
forward conceptually, since it does not require the use of any 
tetrads or proper area calculations. One simply shoots rays back- 
ward from a distant observer and integrates the geodesic path 
until the ray crosses the midplane, where the fluid four-velocity 
can be determined analytically as in Novikov & Thorne (1973). 
This gives the redshift of the emission line as seen by the ob- 
server, and the spectrum is given by the invariant /„/v 3 . 



Figure 7. Error estimate £ as a function of photon number for a relativistically 
broadened iron line, defined in Equation (55). Forty inclination bins were used. 


In Figure 6, we show the shape of a relativistically broadened 
emission line as viewed by observers at different inclination 
angles for the spin values a/M = 0 and 0.99. In all cases, 
the emissivity profile scales like I ~ r~ 3 and the outer disk 
is truncated at r = 15 M. The disk extends all the way into 
the horizon, with the fluid velocity inside the inner-most stable 
circular orbit (ISCO) determined by conserving the energy and 
angular momentum at the ISCO and solving for u r from the 
relation = — 1. For the observer-to-emitter calculation, 

we use the same ray-tracing code described in Schnittman & 
Bertschinger (2004), with 10 7 photons evenly spaced in the 
image plane for each inclination. We find excellent agreement 
in all cases, validating our emitter-to-observer techniques, at 
least for planar test-particle orbits. 

This test, in turn, naturally leads to a simple convergence test 
for our Monte Carlo code. Integrating over energy and observer 
inclination angle i, we can apply the following metric to estimate 
the error due to the use of a finite number of photons: 

_ [f d cos / / dE(I\ 0 (E) — 7 h j(E)) 2 ] 1/2 
[f cl cos i f dEI*(E)] 1/2 

where I\ 0 (E ) is the spectrum calculated at low resolution, com- 
pared with the theoretically perfect spectrum Ihi(E) calculated at 
high resolution. As expected for a Monte Carlo calculation, we 
find that the total error scales with photon number like A -1 / 2 , as 
shown in Figure 7. This is consistent with similar spectral cal- 
culations done with the Monte Carlo radiation code grmonty 
(Dolence et al. 2009). Also shown in Figure 7 are the errors 
e(/V) for the observer-to-emitter approach, using a total of 40 
inclinations for both cases. We note that the emitter-to-observer 
method is more than a factor of two more efficient for the same 
calculation. This is because we can selectively shoot more pho- 
tons from the inner regions of the disk, but in the reverse method, 
the photons are distributed evenly in the image plane (this uni- 
form distribution is not strictly necessary; e.g.. Noble et al. 2007 
use an adaptive grid to improve resolution in bothros). 

The next test is similar, but also includes polarization effects. 
Instead of an emission line with I(r) ~ r -3 , we use the diluted 
thermal spectrum for a Novikov-Thorne (NT) disk with an inner 
edge at the ISCO. The emission has the polarization and limb 
darkening appropriate for a scattering-dominated atmosphere 
(Chandrasekhar 1960). For the observer-to-emitter approach, 
in addition to utilizing the /,,/y 3 invariance, we also parallel 
transport two polarization basis vectors corresponding to the 
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Figure 8. Comparison of observer-to-emitter and emitter-to-observer polarized images for a NT disk with polarization given by a scattering-dominated atmosphere. 
The disk extends from an inner edge at the ISCO out to R 0 ut = 15 M. The black hole has spin a/M = 0.99 and the observer is at an inclination of i = 75°. The 
intensity color scale is logarithmic and the polarization vectors are linearly proportional to the local degree of polarization, as observed at infinity. 

(A color version of this figure is available in the online journal.) 




Figure 9. Polarization degree and angle as a function of photon energy for a black hole with spin a/M = 0.99, luminosity 0.1 L Edd> and mass 10 Mq. The disk extends 
from the ISCO out to r = 15M. The observer-to-emitter (solid curves) and emitter-to-observer (diamonds) methods agree closely over a range of inclinations. 

(A color version of this figure is available in the online journal.) 



Figure 10. Outgoing intensity and degree of polarization for radiation emitted 
corona optical depths. The face-on orientation is unpolarized due to symmetry. 

two axes in the observer plane normal to the photon propagation 
direction. Then, when the ray intersects the disk, we calculate 
a local tetrad in order to determine the local angle of incidence 
and thus the degree of polarization. The direction of polarization 
is projected onto the parallel-transported basis vectors to give 
the observed angle at infinity. 

The two approaches give identical results, for a Kerr black 
hole with spin a/M = 0.99, R oat = 15 M, and observer 
inclination angle 75°, as shown in the images in Figure 8. 
The color code is logarithmic in total intensity and covers four 
orders of magnitude and the small vectors scale linearly with 
the degree of polarization. For the purposes of comparison, we 



0.0 0.2 0.4 0.6 0.8 1.0 

cos i 


a scattering-dominated atmosphere, as a function of inclination, for a range of 

have not included returning radiation here, despite the important 
effect it has on the polarization signal (Agol & Krolik 2000; 
Schnittman & Krolik 2009). In fact, it is precisely due to the 
critical importance of returning radiation that we were forced 
to employ the emitter-to-observer approach in Schnittman & 
Krolik (2009, 2010). 

In Figure 9, we show the observables of polarization degree 
and angle as a function of energy for a range of inclination 
angles, assuming an Eddington-scaled accretion rate of m = 0.1 
and black hole mass M = 10 Mq. Again, we find excellent 
agreement between the emitter-to-observer and observer-to- 
emitter methods. 
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Figure 11. Stokes parameters for scattering of unpolarized light off an optically thick atmosphere. See the text for a description. Compare with Figures 24 and 25 of 
Chandrasekhar (1960). 
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Figure 12. Spectra and polarization of flux from a disk and corona with a slab geometry corresponding to an AGN accretion disk, as described in the text. The solid 
curves correspond to the total flux, while the dotted, dashed, dot-dashed, triple-dot-dashed, and long-dashed curves correspond to A scat = 0, 1, 2, 3, and ^ 5. For 
clarity, only the solid curves are shown for the polarization degree. Compare with Figure 5 of Poutanen & Svensson (1996). 



Next, we move on to testing the coronal scattering algorithms. 
We focus on a plane-parallel geometry with an optically thick 
disk covered by a corona with variable optical depth r and 
electron temperature T e . In Figure 10, we show the effects of a 
scattering atmosphere on the emergent flux and polarization as 
a function of angle. The seed photons are emitted isotropically 
(i.e., with a cos 0 weight appropriate for an optically thick 
disk) from the disk surface with zero polarization, then scatter 
through a cold corona. Photons that scatter back to the disk are 
reflected via the diffuse scattering formula of Equation (53). 
In the limit of r —> oo, we reproduce the limb darkening 
and horizontal polarization results from Chandrasekhar (1960), 
Table XXIV. 

In Figure 1 1 , we carry out a similar scattering experiment, 
but with the seed photons incident from above the disk along 


a single direction. Setting r = 10, we should reproduce the 
analytic diffuse reflection expressions of Chandrasekhar (1960), 
shown in that paper in Figures 24 and 25 for an incident 
unpolarized beam with cos 0o = 0.8, 0.5, and <fo = 0. Following 
Chandrasekhar (1960), we plot the Stokes parameters /, Q, and 
U as a function of reflection angle, normalized to the incident 
intensity. The asterisks are the Monte Carlo calculations and the 
solid curves are the analytic predictions. 

On the left-hand side of each plot, we show the polarization 
as a function of 6 for ip — (po = 0°, ±180°. The value of 9q 
is designated with a vertical dashed line. Negative values of 0 
correspond to photons reflected back in the general direction 
of the incident photons, i.e., (p — cpo = 0. Thus, we see a 
natural peak in the intensity corresponding to backscattering 
as in Equation (42). Similarly, the degree of polarization is 
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Figure 13. Same as Figure 12, but for T e = 352 keV and r = 0.05. Compare with Figure 6 of Poutanen & Svensson (1996). 


maximized for 90° scattering and oriented in the plane of the 

disk (Q > 0). 

On the right-hand side of each plot, we show the Stokes 
parameters for photons scattered with (p — (pa = ±90°. In this 
case, the planar symmetry is broken and we find a non-zero 
value of U. Again, the degree of polarization is maximized for 
scattering angles near 90°. 

Lastly, we test the inverse Compton effects of a hot corona 
by reproducing the AGN-type spectra of Poutanen & Svensson 
(1996). The seed photons are again isotropic and unpolarized, 
with a blackbody spectrum with 7bb = 10 eV. When reflecting 
off the cold disk, we implement the Compton recoil losses of 
Equation (54). Following Poutanen & Svensson (1996), we also 
include atomic absorption in the disk with an extremely simple 
toy model based on the photoelectric cross sections of Morrison 
& McCammon (1983). 

In Figure 1 2, the corona temperature is 56 keV with an optical 
depth r = 0.5. In Figure 13, T e = 352 keV and r = 0.05, 
corresponding to Figures 5 and 6 in Poutanen & Svensson 
(1996). In the top panels, we show the observed flux at two 
inclination angles cos i = 0.11,0.5 and in the bottom panels 
we show the polarization degree S(% ) = Q/I x 100. In all 
panels, the solid curves correspond to the total flux, while the 
dotted, dashed, dot-dashed, triple-dot-dashed, and long-dashed 
curves represent subsets of the flux, binned by the number of 
coronal scatterings 0, 1,2, 3, and Jj 5. Photon packets that return 
to the disk suffer photoelectric absorption and Compton recoil 
losses and are then launched again from the disk, resetting N scat 
to zero. Thus, the dotted curves in Figures 12 and 13 have 
significant power around the Compton hump at 10-100 keV. 
As discussed in Schnittman & Krolik (2010), more scatterings 
in a sandwich corona effectively constrain the geometry and 
increase the amplitude of polarization at high energies. 


We find excellent agreement overall, but the spectra are 
clearly dominated by Monte Carlo noise above ~ 1 00 keV. For 
these disk and coronal parameters, this corresponds to seed 
photons that have already scattered on average over 25 times, 
so it is very difficult to resolve any polarization signal at the few 
percent level. Additionally, due to our photon packet algorithm, 
we are limited to energy-independent electron cross sections, so 
we should expect that the accuracy of our spectral predictions 
breaks down much above 100 keV anyway. 

6. CONCLUSIONS 

We have presented the technical details behind the general 
relativistic radiation transport code Pandurata. Its capabilities 
include optically thin emission and absorption, Compton scat- 
tering, polarization, spectral and timing analysis, and flexible 
geometries that allow analysis of numerous accretion models 
and MHD simulations. We have discussed a number of practi- 
cal challenges that other teams may also face when working to 
develop similar ray-tracing codes, such as the method of weights 
in the scattering kernel. 

This is by no means the final word on Pandurata. Its great 
strength lies in its flexibility and we envision numerous up- 
grades and improvements in the near future. These will include, 
but not be limited to, detailed ionization balance in the disk 
photosphere for improved AGN modeling. Improved spectra 
accuracy for atomic absorption and emission will likely require 
a more traditional mono-energetic photon packet architecture, 
as done in Dolence et al. (2009). This should be quite straight- 
forward to incorporate into Pandurata, and will also allow us to 
include Klein-Nishina cross sections at high energies, albeit at 
the loss of computational efficiency that comes with broadband 
photon packets. 
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Time interpolation between simulation snapshots for generat- 
ing more accurate light curves will naturally allow for relativistic 
effects such as Lorentz contraction and time delays for emission 
close to the black hole (Noble & Krolik 2009). Future timing 
applications would look at soft/hard X-ray lags from MHD 
simulations, energy-dependent, quasi-periodic oscillations, and 
broadband noise. The inclusion of more sophisticated emission 
and absorption processes (e.g., angle-dependent synchrotron) 
will allow us to model low-luminosity and radiatively inefficient 
sources such as Sgr A* and M87. Perhaps most importantly, we 
will work to close the final remaining gap between theory and 
observation by incorporating Pandurata spectra into a data 
analysis framework like xspec and making it publicly available 
to the X-ray astronomy community. 


The momentum equations are a bit more involved, but there are 
only two of them (for p r and po', p,j, is conserved): 


dp r 

dt 


dco p, + wpj, da 
Pip + — + 


dr 


a dr 2(p,+wp<p) 


A , 1 , 1 , 

o I —Pr + —Pe + —Pip 

' p z p z U7 Z 


8 

dr 


(A5a) 


dp e 
dt 


dco p t + copA, da 
~P<P + — + 


30 ■ 


30 


2{p, + cop 0 ) 

1 


7T77 — Pr + —P~6 + —P* 

30 \p z p z TU~ v 


(A5b) 
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APPENDIX A 

HAMILTONIAN EQUATIONS OF MOTION 

The equations of motion for the Hamiltonian H in 
Boyer-Lindquist coordinates, as given in Section 3, are repeated 
here for completeness: 


H(r , 9, cp, p r , pg, p</,) 



= -Pt = coprf, + a ( 

( ^ 2 1 2 1 2 2^ 
^P> + S Pd + ^ P * +m p 

r 

and according to classical theory: 



dx l dH 

dt dpi 

(Ala) 


dpt dH 

dt dx‘ 

(Alb) 


For convenience of notation, we define the quantity D 2 as 


D 2 (r , 0, cp, p r , pg, p<j,) = -p 2 + ^-pg + — -pl+m 2 . (A2) 

p z p z m L r 

Then, for an arbitrary variable y e (x l , pi ), the partial derivative 
of H can be written 


dH d da 

^ = 87 ( "'’‘ ,+ s7 D 


dD 2 


2 p, + cop $ dy 


(A3) 


The first set of Hamiltonian’s equations are straightforward tc 
produce: 


The relevant spatial derivatives are as follows: 


dco co 2 


dr 2 Ma 


3 r 2 + a 2 { 1 + cos 2 0) cos 2 6 

r z 


dco co 


2 r 


dO 2 Ma 


y2Ma 2 — a 2 r — — ) sin 0 cos 0 


da 1 da 2 

dr 2 a dr 

da 1 da 2 

90 ~ 2a~W ’ 


da 2 
dr 

da 2 

~dQ 


2 M\(a A -r A 2r 2 a 2 sin 2 0 


A p 2 ) \ A p z 

AMci 2 r sin0 cos0(a 2 + r 2 ) 


d ( 1 


dr \uj 2 


sin 2 6 1 r + 


A p 2 

2 Ma 2 sin 2 0(a 2 cos 2 6 


(A6a) 

(A6b) 

(A6c) 

(A6d) 

, (A6e) 

(A6f) 
■ r 2 )^" 1 


(A6g) 


-f 1 

dO l m 2 


4 sin 0 cos 6 




+ (r 2 + a 2 ) 


2Mcr sin 2 9 


r 2 + a 2 1 

p 4 + p 2 




dr \p 


(A6h) 


(A6i) 


dr 

dH x 

p r a 2 A 

(A4a) 

3 

dt 

dp r 

p, + cop 0 p 2 

3 9 

d9 

3//, 

Pe a 2 

(A4b) 

3 

dt 

3 Pg 

Pt + u>P<p P 2 ’ 

dr 

dtp 

dH\ 

P0 a 2 

(A4c) 

3 

dt 

dp<P 

p, + copt/f m 2 

3 9 



= —a~ A sin 6 cos 6 , 
P 4 


2r 


= —a 1 sin 9 cos 9 . 
P 4 


(A6j) 

(A6k) 
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APPENDIX B 

MONTE CARLO SAMPLING OF THE 
MAXWELL-JUTTNER DISTRIBUTION 

For any normalized distribution function f(x ) with x e 
(— oo, oo), one can always define the cumulative distribution 
function 

cdf(x) = F(x) = / f(x')dx', (Bl) 

J — OO 

with F(—o o) = 0 and F(o o) = 1. Then, by selecting a uniform 
random number k e [0, 1), the choice x = F '( a ) will be 
distributed according to fix). However, in most cases, Fix) 
cannot be written in closed form, so other methods are required. 

One simple technique described in Press et al. (1992) is the 
“rejection method,” where an auxiliary function g(x) is used, 
where g(x) > f(x) everywhere and G(x) is easy to calculate. 
We begin by selecting a trial xo = G ~ 1 (/-o), then pick another 
random deviate M. If M < f{xf)/g(xf), thenxo is selected as 
a representative sample of /(x); otherwise, we try again with a 
new A. 0 . Of course, if g(x) is large enough, it is easy to ensure 
that it is greater than f(x) everywhere. However, the efficiency 
of this method is limited by the ratio of the areas under the two 
curves f(x) and g(x), so it is desirable to pick g(x) as close to 
fix) as possible (Press et al. 1992). 

For the Maxwell-Iuttner distribution defined in Equa- 
tion (38): 

/(y)~ y 2 /3exp(-y/0 r ) (B2) 

we choose an auxiliary function 


For angles uniformly distributed in cos 9 and (f>, one can show 
that the probability distribution function (pdf) for w is 


P(w)= ^ 



(C3) 


for 0 ^ w ^ 3/2 and 0 otherwise. 

For multiply-scattered photons, the weight function is multi- 
plicative, since the individual scattering events are uncorrelated. 
For n scatters, the net weight is given by 

n 

W = Y\wi. (C4) 

1 = 1 

To determine the pdf P(W ), we define a new variable Z: 

n n 

Z = In W = ^ In wi = ^ Zi . (C5) 

! = 1 1=1 


For large values of n, the central limit theorem dictates that the 
distribution of Z should be Gaussian: 


P(Z) = 


1 / (Z-np, z ) 2 \ 

a z jinn ^ \ 2 no} / ’ 


(C6) 


where /i z and err are the mean and variance of P(z), respectively. 
From Equation (C3) and the variable transformation z = In w, 
we have 


g(y) ~ y 2 exp(-y/6 T ). 


This gives 


G(y) = 1 - 


e - y / 0 T 29 2 +26 T y + y 2 
e~ l 'G 2 9 2 +29 t + \ 


(B3) 


(B4) 


P( Z ) = P ( w ) 


dw 

dz 



(C7) 


with z e (—oo, In 3/2]. This gives ii z = —0.208 and a 2 = 
0.710. Now, we see that the pdf P(W) is given by a log-normal 
distribution: 


for the cumulative distribution function. Inverting (B4) is not 
trivial, but can be done numerically with a simple root finder. For 
these choices of f(y) and g(y), we find an excellent efficiency 
for this algorithm of ~90%. 

APPENDIX C 

COMPARISON OF SCATTERING KERNELS 


1 

P(W) = -j= exp 

Wa z ^/2jtn 


(In W — nfx z ) 2 \ 

2 no} J 


(C8) 


For photons random walking through an atmosphere of 
optical depth r, from numerical experiments we find the pdf 
of the number of scatters required to escape can be closely 
approximated by 


As described in Section 4, there are (at least) two different 
ways to implement the scattering of polarized light off of free 
electrons. 

The method of weights picks a random scattering angle from 
a uniform distribution of cos 0 e [—1, 1] and </> e [0, 2 n), then 
weights the scattered beam of photons by the cross section in that 
direction, normalized by the average cross section to conserve 
flux. By integrating Equation (43) over <p , this resembles the 
classical cross section for unpolarized light: 

/' 3 

w(9) = — = -(cos 2 9 + 1). (Cl) 

I 4 

Because repeated scatters tend to increase the level of po- 
larization (indeed, in the microscopic limit, every photon has 
S — 1 ), we will focus on the case where <5=1, giving 

3 

w(9, (p) — -(cos 2 9 cos 2 f + sin 2 </>). (C2) 


n / n \ 

P{n)= ^ eXp (“^)' (C9) 
Then, the net distribution P(W) for all scatting orders is simply 

POO 

P{W\x)= \ dn Pin; r) P(W\ n). (CIO) 
Jo 

The relative contribution to the spectrum from photons with a 
weight in the range ( W , W + dW ) is P(W)WdW , so we require 
PfW) to decrease faster than W ~ 2 for large W if the calculation 
is to converge. In Figure 14, we plot W 2 P(W) for a range of r. 
Our analytic results suggest that for r > 2, any polarization 
spectrum formed using this Monte Carlo weighting method 
should be dominated by the rarest, highest-weight photon 
packets, confirming what we have seen in trial runs with large r . 
Now, in practice, the convergence is not quite as bad as Figure 14 
suggests, for two primary reasons. First, the seed photon 
packets have little or no polarization, so the initial weighting 
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w 


Figure 14. Relative contribution to the total observed spectrum by photons of a 
given weight, for optical depths of r = 1, 2, 3, 5, and 10. Any calculation with 
r > 2 will not formally converge using the method of weights. 

function more closely resembles Equation (Cl), which leads 
to a significantly tighter range in W: /x unpo i = —0.027 and 
^unpoi = 0.047 (in this unpolarized limit, the weight method 
converges for all optical depths up to r > 200). Second, for the 
small-to-moderate optical depths of r < 5, the typical number 
of scatters n is still small enough that the mean value theorem 
does not strictly apply, in effect cutting off the high-weight tails 
in Equation (C6) and further reducing the contribution from 
statistical outliers. 

However, for r > 10, the polarization of a typical photon 
bundle reaches 8 -> 1 after just a few scatters and the large 
number of total scattering events allows us to reproduce these 
analytic results with numerical tests of the Monte Carlo code. In 
Figure 15, we show the distribution of weights from a calculation 
using unpolarized seed photons, scattering through optical 
depths of r = 2, 5, and 10. While we find that the r = 5 case 
does converge eventually, in practice we find the convergence 
is so slow that another Monte Carlo method is preferable. 
Furthermore, the highest weights have the fewest events and 
thus also suffer from small-number statistics, potentially adding 
to the “undue influence” of outliers. This can be seen in the 
scatter at the high-weight end of each dataset. 

Instead of picking a scattering angle at random and weighting 
it by Equation (47), let us use the differential cross section 
(Equation (43)) to obtain the scattering pdf: 

P(0,0)=-^-[(l-5)(cos 2 0 + l) 

l07T 

+ 28 cos 2 0 cos 2 0 + 28 sin 2 </>]. (Cll) 

Integrating over 0, we again find the standard Thomson cross 
section, which holds even for <5^0: 

3 9 

P(cos<9) = -(cos 2 (9 + 1). (Cl 2) 

8 

Writing /x = cos 6 for convenience, the cumulative distribution 
function is given by 

cdf( M ) = R(/xW = + 3/x + 4) (Cl 3) 

To pick an appropriate value for /x, we generate a random 
number X from a uniform distribution X e [0, 1) and invert 



W 


Figure 15. Relative contribution to the total observed spectrum by photons of a 
given weight, for optical depths of r = 2, 5, and 10, sorted by color (red, blue, 
and black). The solid lines are the analytic results and the crosses are “data” 
from Monte Carlo calculations of 10 6 photons each. 

(A color version of this figure is available in the online journal.) 

Equation (Cl 3), in effect solving for the root of the cubic: 

/X 3 + 3/x + 4 - 8k = 0. (C14) 

Because cdf(/x) is monotonically increasing, this equation is 
guaranteed to have a single real root in the interval — 1 ^ /x ^ 1 . 

Once /x is selected, we choose 0 by the same method, now 
using the pdf 

P(0; /x) = \ [(1 - 5)(/x 2 + 1) 

27T(fl l + 1 ) 

+ 28 /J? cos 2 0 + 2<$ sin 2 0] , (C15) 

which gives 

0 8 ( u? — 1 \ 

cdf(0; /x) = — sin(2 0). (C16) 

2n Ait \ /x- + 1 / 

Again, to pick an appropriate 0 given a uniform random X, one 
must invert Equation (Cl 6) to obtain 0 = cdf~ l (7,j. Unfortu- 
nately, this is equivalent to solving Kepler’s equation, which has 
no closed-form solution and must be done numerically. Fortu- 
nately, this is equivalent to solving Kepler’s equation, one of the 
best-studied numerical problems in astrophysics! In practice, 
we use the iterative approach outlined in Murray & Dermott 
(1999). While slightly more time consuming than the method of 
weights, the exact cross section method has the distinct advan- 
tage of converging for an arbitrary number of scatterings and 
thus is the method we prefer for Pandurata. 
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